home *** CD-ROM | disk | FTP | other *** search
/ Aminet 30 / Aminet 30 (1999)(Schatztruhe)[!][Apr 1999].iso / Aminet / gfx / misc / gnuplot-3.7src.lha / gnuplot-3.7src / gnuplot-3.7.lha / gnuplot-3.7 / interpol.c < prev    next >
C/C++ Source or Header  |  1998-11-16  |  30KB  |  1,031 lines
  1. #ifndef lint
  2. static char *RCSid = "$Id: interpol.c,v 1.29 1998/04/14 00:15:45 drd Exp $";
  3. #endif
  4.  
  5. /* GNUPLOT - interpol.c */
  6.  
  7. /*[
  8.  * Copyright 1986 - 1993, 1998   Thomas Williams, Colin Kelley
  9.  *
  10.  * Permission to use, copy, and distribute this software and its
  11.  * documentation for any purpose with or without fee is hereby granted,
  12.  * provided that the above copyright notice appear in all copies and
  13.  * that both that copyright notice and this permission notice appear
  14.  * in supporting documentation.
  15.  *
  16.  * Permission to modify the software is granted, but not the right to
  17.  * distribute the complete modified source code.  Modifications are to
  18.  * be distributed as patches to the released version.  Permission to
  19.  * distribute binaries produced by compiling modified sources is granted,
  20.  * provided you
  21.  *   1. distribute the corresponding source modifications from the
  22.  *    released version in the form of a patch file along with the binaries,
  23.  *   2. add special version identification to distinguish your version
  24.  *    in addition to the base release version number,
  25.  *   3. provide your name and address as the primary contact for the
  26.  *    support of your modified version, and
  27.  *   4. retain our contact information in regard to use of the base
  28.  *    software.
  29.  * Permission to distribute the released version of the source code along
  30.  * with corresponding source modifications in the form of a patch file is
  31.  * granted with same provisions 2 through 4 for binary distributions.
  32.  *
  33.  * This software is provided "as is" without express or implied warranty
  34.  * to the extent permitted by applicable law.
  35. ]*/
  36.  
  37.  
  38. /*
  39.  * C-Source file identification Header
  40.  *
  41.  * This file belongs to a project which is:
  42.  *
  43.  * done 1993 by MGR-Software, Asgard  (Lars Hanke)
  44.  * written by Lars Hanke
  45.  *
  46.  * Contact me via:
  47.  *
  48.  *  InterNet: mgr@asgard.bo.open.de
  49.  *      FIDO: Lars Hanke @ 2:243/4802.22   (as long as they keep addresses)
  50.  *
  51.  **************************************************************************
  52.  *
  53.  *   Project: gnuplot
  54.  *    Module:
  55.  *      File: interpol.c
  56.  *
  57.  *   Revisor: Lars Hanke
  58.  *   Revised: 26/09/93
  59.  *  Revision: 1.0
  60.  *
  61.  **************************************************************************
  62.  *
  63.  * LEGAL
  64.  *  This module is part of gnuplot and distributed under whatever terms
  65.  *  gnuplot is or will be published, unless exclusive rights are claimed.
  66.  *
  67.  * DESCRIPTION
  68.  *  Supplies 2-D data interpolation and approximation routines
  69.  *
  70.  * IMPORTS
  71.  *  plot.h
  72.  *    - cp_extend()
  73.  *    - structs: curve_points, coordval, coordinate
  74.  *
  75.  *  setshow.h
  76.  *    - samples, is_log_x, base_log_x, xmin, xmax, autoscale_lx
  77.  *    - plottypes
  78.  *
  79.  *  proto.h
  80.  *    - solve_tri_diag()
  81.  *    - typedef tri_diag
  82.  *
  83.  * EXPORTS
  84.  *  gen_interp()
  85.  *  sort_points()
  86.  *  cp_implode()
  87.  *
  88.  * BUGS and TODO
  89.  *  I would really have liked to use Gershon Elbers contouring code for
  90.  *  all the stuff done here, but I failed. So I used my own code.
  91.  *  If somebody is able to consolidate Gershon's code for this purpose
  92.  *  a lot of gnuplot users would be very happy - due to memory problems.
  93.  *
  94.  **************************************************************************
  95.  *
  96.  * HISTORY
  97.  * Changes:
  98.  *  Nov 24, 1995  Markus Schuh (M.Schuh@meteo.uni-koeln.de):
  99.  *      changed the algorithm for csplines
  100.  *      added algorithm for approximation csplines
  101.  *      copied point storage and range fix from plot2d.c
  102.  *
  103.  *  Dec 12, 1995 David Denholm
  104.  *      oops - at the time this is called, stored co-ords are
  105.  *      internal (ie maybe log of data) but min/max are in
  106.  *      user co-ordinates. 
  107.  *      Work with min and max of internal co-ords, and
  108.  *      check at the end whether external min and max need to
  109.  *      be increased. (since samples is typically 100 ; we
  110.  *      dont want to take more logs than necessary)
  111.  *      Also, need to take into account which axes are active
  112.  *
  113.  *  Jun 30, 1996 Jens Emmerich
  114.  *      implemented handling of UNDEFINED points
  115.  */
  116.  
  117. #include "plot.h"
  118. #include "setshow.h"
  119.  
  120.  
  121. /* in order to support multiple axes, and to simplify ranging in
  122.  * parametric plots, we use arrays to store some things. For 2d plots,
  123.  * elements are z=0,y1=1,x1=2,z2=4,y2=5,x2=6 these are given symbolic
  124.  * names in plot.h
  125.  */
  126.  
  127. extern double min_array[AXIS_ARRAY_SIZE];
  128. extern double max_array[AXIS_ARRAY_SIZE];
  129. extern int auto_array[AXIS_ARRAY_SIZE];
  130. extern TBOOLEAN log_array[AXIS_ARRAY_SIZE];
  131. extern double base_array[AXIS_ARRAY_SIZE];
  132. extern double log_base_array[AXIS_ARRAY_SIZE];
  133.  
  134.  
  135. #define Inc_c_token if (++c_token >= num_tokens)    \
  136. int_error ("Syntax error", c_token);
  137.  
  138.  
  139. /*
  140.  * IMHO, code is getting too cluttered with repeated chunks of
  141.  * code. Some macros to simplify, I hope.
  142.  *
  143.  * do { } while(0) is comp.lang.c recommendation for complex macros
  144.  * also means that break can be specified as an action, and it will
  145.  * 
  146.  */
  147.  
  148.  
  149. /* store VALUE or log(VALUE) in STORE, set TYPE as appropriate
  150.  * Do OUT_ACTION or UNDEF_ACTION as appropriate
  151.  * adjust range provided type is INRANGE (ie dont adjust y if x is outrange
  152.  * VALUE must not be same as STORE
  153.  */
  154.  
  155. #define STORE_AND_FIXUP_RANGE(STORE, VALUE, TYPE, MIN, MAX, AUTO, OUT_ACTION, UNDEF_ACTION)\
  156. do { STORE=VALUE; \
  157.     if (TYPE != INRANGE) break;  /* dont set y range if x is outrange, for example */ \
  158.     if ( VALUE<MIN ) { \
  159.        if (AUTO & 1) MIN=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; break; }  \
  160.     } \
  161.     if ( VALUE>MAX ) {\
  162.        if (AUTO & 2) MAX=VALUE; else { TYPE=OUTRANGE; OUT_ACTION; }   \
  163.     } \
  164. } while(0)
  165.  
  166. #define UPDATE_RANGE(TEST,OLD,NEW,AXIS) \
  167. do { if (TEST) { \
  168.      if (log_array[AXIS]) OLD = pow(base_array[AXIS], NEW); else OLD = NEW; \
  169.      } \
  170. } while(0)
  171.  
  172. /* use this instead empty macro arguments to work around NeXT cpp bug */
  173. /* if this fails on any system, we might use ((void)0) */
  174.  
  175. #define NOOP            /* */
  176.  
  177. #define spline_coeff_size 4
  178. typedef double spline_coeff[spline_coeff_size];
  179. typedef double five_diag[5];
  180.  
  181. static int next_curve __PROTO((struct curve_points * plot, int *curve_start));
  182. static int num_curves __PROTO((struct curve_points * plot));
  183. static double *cp_binomial __PROTO((int points));
  184. GP_INLINE static double s_pow __PROTO((double base, unsigned int exponent));
  185. static void eval_bezier __PROTO((struct curve_points * cp, int first_point, int num_points, double sr, coordval * px, coordval * py, double *c));
  186. static void do_bezier __PROTO((struct curve_points * cp, double *bc, int first_point, int num_points, struct coordinate * dest));
  187. static int solve_five_diag __PROTO((five_diag m[], double r[], double x[], int n));
  188. static spline_coeff *cp_approx_spline __PROTO((struct curve_points * plot, int first_point, int num_points));
  189. static spline_coeff *cp_tridiag __PROTO((struct curve_points * plot, int first_point, int num_points));
  190. static void do_cubic __PROTO((struct curve_points * plot, spline_coeff * sc, int first_point, int num_points, struct coordinate * dest));
  191. static int compare_points __PROTO((struct coordinate * p1, struct coordinate * p2));
  192.  
  193.  
  194. /*
  195.  * position curve_start to index the next non-UNDEFINDED point,
  196.  * start search at initial curve_start,
  197.  * return number of non-UNDEFINDED points from there on,
  198.  * if no more valid points are found, curve_start is set
  199.  * to plot->p_count and 0 is returned
  200.  */
  201.  
  202. static int next_curve(plot, curve_start)
  203. struct curve_points *plot;
  204. int *curve_start;
  205. {
  206.     int curve_length;
  207.  
  208.     /* Skip undefined points */
  209.     while (*curve_start < plot->p_count
  210.        && plot->points[*curve_start].type == UNDEFINED) {
  211.     (*curve_start)++;
  212.     };
  213.     curve_length = 0;
  214.     /* curve_length is first used as an offset, then the correkt # points */
  215.     while ((*curve_start) + curve_length < plot->p_count
  216.       && plot->points[(*curve_start) + curve_length].type != UNDEFINED) {
  217.     curve_length++;
  218.     };
  219.     return (curve_length);
  220. }
  221.  
  222.  
  223. /* 
  224.  * determine the number of curves in plot->points, separated by
  225.  * UNDEFINED points
  226.  */
  227.  
  228. static int num_curves(plot)
  229. struct curve_points *plot;
  230. {
  231.     int curves;
  232.     int first_point;
  233.     int num_points;
  234.  
  235.     first_point = 0;
  236.     curves = 0;
  237.     while ((num_points = next_curve(plot, &first_point)) > 0) {
  238.     curves++;
  239.     first_point += num_points;
  240.     }
  241.     return (curves);
  242. }
  243.  
  244.  
  245.  
  246. /*
  247.  * build up a cntr_struct list from curve_points
  248.  * this funtion is only used for the alternate entry point to
  249.  * Gershon's code and thus commented out
  250.  ***deleted***
  251.  */
  252.  
  253.  
  254. /*
  255.  * cp_binomial() computes the binomial coefficients needed for BEZIER stuff
  256.  *   and stores them into an array which is hooked to sdat.
  257.  * (MGR 1992)
  258.  */
  259. static double *cp_binomial(points)
  260. int points;
  261. {
  262.     register double *coeff;
  263.     register int n, k;
  264.     int e;
  265.  
  266.     e = points;            /* well we're going from k=0 to k=p_count-1 */
  267.     coeff = (double *) gp_alloc(e * sizeof(double), "bezier coefficients");
  268.  
  269.     n = points - 1;
  270.     e = n / 2;
  271.     coeff[0] = 1.0;
  272.  
  273.     for (k = 0; k < e; k++) {
  274.     coeff[k + 1] = coeff[k] * ((double) (n - k)) / ((double) (k + 1));
  275.     }
  276.  
  277.     for (k = n; k >= e; k--)
  278.     coeff[k] = coeff[n - k];
  279.  
  280.     return (coeff);
  281. }
  282.  
  283.  
  284. /* This is a subfunction of do_bezier() for BEZIER style computations.
  285.  * It is passed the stepration (STEP/MAXSTEPS) and the addresses of
  286.  * the double values holding the next x and y coordinates.
  287.  * (MGR 1992)
  288.  */
  289.  
  290. /*
  291.  * well, this routine runs faster with the 68040 striptease FPU
  292.  * and it handles zeroes correctly - there had been some trouble with TC
  293.  */
  294.  
  295. GP_INLINE static double s_pow(base, exponent)
  296. double base;
  297. unsigned int exponent;
  298. {
  299.     double y;
  300.  
  301.     if (exponent == 0)
  302.     return (1.0);
  303.     if (base == 0.0)
  304.     return (0.0);
  305.  
  306.     /* consider i in binary = abcd
  307.      * x^i = x^(8a+4b+2c+d) = x^8a * x^4b * x^2b * x^d
  308.      *                      = a?x^2^2^2:1 * b?x^2^2:1 + ...
  309.      * so for each bit in exponent, square x, multiplying it into accumulator
  310.      *
  311.      */
  312.  
  313.     y = 1;
  314.     while (exponent) {
  315.     if (exponent & 1)
  316.         y *= base;
  317.     base *= base;
  318.     /* if exponent was signed, this could be trouble ! */
  319.     exponent >>= 1;
  320.     }
  321.  
  322.     return (y);
  323. }
  324.  
  325. static void eval_bezier(cp, first_point, num_points, sr, px, py, c)
  326. struct curve_points *cp;
  327. int first_point;        /* where to start in plot->points (to find x-range) */
  328. int num_points;            /* to determine end in plot->points */
  329. double sr;
  330. coordval *px;
  331. coordval *py;
  332. double *c;
  333. {
  334.     unsigned int n = num_points - 1;
  335.     /* HBB 980308: added 'GPHUGE' tag for DOS */
  336.     struct coordinate GPHUGE *this_points;
  337.  
  338.     this_points = (cp->points) + first_point;
  339.  
  340.     if (sr == 0.0) {
  341.     *px = this_points[0].x;
  342.     *py = this_points[0].y;
  343.     } else if (sr == 1.0) {
  344.     *px = this_points[n].x;
  345.     *py = this_points[n].y;
  346.     } else {
  347.     unsigned int i;
  348.     double lx = 0.0, ly = 0.0;
  349.     double sr_to_the_i = 1.0;
  350.     double dsr_to_the_di = s_pow(1 - sr, n);
  351.     double reciprocal_dsr = 1.0 / (1 - sr);
  352.  
  353.     for (i = 0; i <= n; i++) {
  354.         double u = c[i] * sr_to_the_i * s_pow(1 - sr, n - i);
  355.         lx += this_points[i].x * u;
  356.         ly += this_points[i].y * u;
  357.         sr_to_the_i *= sr;
  358.         dsr_to_the_di *= reciprocal_dsr;
  359.     }
  360.  
  361.     *px = lx;
  362.     *py = ly;
  363.     }
  364. }
  365.  
  366. /*
  367.  * generate a new set of coordinates representing the bezier curve and
  368.  * set it to the plot
  369.  */
  370.  
  371. static void do_bezier(cp, bc, first_point, num_points, dest)
  372. struct curve_points *cp;
  373. double *bc;
  374. int first_point;        /* where to start in plot->points */
  375. int num_points;            /* to determine end in plot->points */
  376. struct coordinate *dest;    /* where to put the interpolated data */
  377. {
  378.     int i;
  379.     coordval x, y;
  380.     int xaxis = cp->x_axis;
  381.     int yaxis = cp->y_axis;
  382.  
  383.     /* min and max in internal (eg logged) co-ordinates. We update
  384.      * these, then update the external extrema in user co-ordinates
  385.      * at the end.
  386.      */
  387.  
  388.     double ixmin, ixmax, iymin, iymax;
  389.     double sxmin, sxmax, symin, symax;    /* starting values of above */
  390.  
  391.     if (log_array[xaxis]) {
  392.     ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  393.     ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  394.     } else {
  395.     ixmin = sxmin = min_array[xaxis];
  396.     ixmax = sxmax = max_array[xaxis];
  397.     }
  398.  
  399.     if (log_array[yaxis]) {
  400.     iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  401.     iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  402.     } else {
  403.     iymin = symin = min_array[yaxis];
  404.     iymax = symax = max_array[yaxis];
  405.     }
  406.  
  407.     for (i = 0; i < samples; i++) {
  408.     eval_bezier(cp, first_point, num_points, (double) i / (double) (samples - 1), &x, &y, bc);
  409.  
  410.     /* now we have to store the points and adjust the ranges */
  411.  
  412.     dest[i].type = INRANGE;
  413.     STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
  414.     STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
  415.  
  416.     dest[i].xlow = dest[i].xhigh = dest[i].x;
  417.     dest[i].ylow = dest[i].yhigh = dest[i].y;
  418.  
  419.     dest[i].z = -1;
  420.     }
  421.  
  422.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  423.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  424.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  425.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  426. }
  427.  
  428. /*
  429.  * call contouring routines -- main entry
  430.  */
  431.  
  432. /* 
  433.  * it should be like this, but it doesn't run. If you find out why,
  434.  * contact me: mgr@asgard.bo.open.de or Lars Hanke 2:243/4802.22@fidonet
  435.  *
  436.  * Well, all this had originally been inside contour.c, so maybe links
  437.  * to functions and of contour.c are broken.
  438.  * ***deleted***
  439.  * end of unused entry point to Gershon's code
  440.  *
  441.  */
  442.  
  443. /* 
  444.  * Solve five diagonal linear system equation. The five diagonal matrix is
  445.  * defined via matrix M, right side is r, and solution X i.e. M * X = R.
  446.  * Size of system given in n. Return TRUE if solution exist.
  447.  *  G. Engeln-Muellges/ F.Reutter: 
  448.  *  "Formelsammlung zur Numerischen Mathematik mit Standard-FORTRAN-Programmen"
  449.  *  ISBN 3-411-01677-9
  450.  *
  451.  * /  m02 m03 m04   0   0   0   0    .       .       . \   /  x0  \    / r0  \
  452.  * I  m11 m12 m13 m14   0   0   0    .       .       . I   I  x1  I   I  r1  I
  453.  * I  m20 m21 m22 m23 m24   0   0    .       .       . I * I  x2  I = I  r2  I
  454.  * I    0 m30 m31 m32 m33 m34   0    .       .       . I   I  x3  I   I  r3  I
  455.  *      .   .   .   .   .   .   .    .       .       .        .        .
  456.  * \                           m(n-3)0 m(n-2)1 m(n-1)2 /   \x(n-1)/   \r(n-1)/
  457.  * 
  458.  */
  459. static int solve_five_diag(m, r, x, n)
  460. five_diag m[];
  461. double r[], x[];
  462. int n;
  463. {
  464.     int i;
  465.     five_diag *hv;
  466.  
  467.     hv = (five_diag *) gp_alloc((n + 1) * sizeof(five_diag), "five_diag help vars");
  468.  
  469.     hv[0][0] = m[0][2];
  470.     if (hv[0][0] == 0) {
  471.     free(hv);
  472.     return FALSE;
  473.     }
  474.     hv[0][1] = m[0][3] / hv[0][0];
  475.     hv[0][2] = m[0][4] / hv[0][0];
  476.  
  477.     hv[1][3] = m[1][1];
  478.     hv[1][0] = m[1][2] - hv[1][3] * hv[0][1];
  479.     if (hv[1][0] == 0) {
  480.     free(hv);
  481.     return FALSE;
  482.     }
  483.     hv[1][1] = (m[1][3] - hv[1][3] * hv[0][2]) / hv[1][0];
  484.     hv[1][2] = m[1][4] / hv[1][0];
  485.  
  486.     for (i = 2; i <= n - 1; i++) {
  487.     hv[i][3] = m[i][1] - m[i][0] * hv[i - 2][1];
  488.     hv[i][0] = m[i][2] - m[i][0] * hv[i - 2][2] - hv[i][3] * hv[i - 1][1];
  489.     if (hv[i][0] == 0) {
  490.         free(hv);
  491.         return FALSE;
  492.     }
  493.     hv[i][1] = (m[i][3] - hv[i][3] * hv[i - 1][2]) / hv[i][0];
  494.     hv[i][2] = m[i][4] / hv[i][0];
  495.     }
  496.  
  497.     hv[0][4] = 0;
  498.     hv[1][4] = r[0] / hv[0][0];
  499.     for (i = 1; i <= n - 1; i++) {
  500.     hv[i + 1][4] = (r[i] - m[i][0] * hv[i - 1][4] - hv[i][3] * hv[i][4]) / hv[i][0];
  501.     }
  502.  
  503.     x[n - 1] = hv[n][4];
  504.     x[n - 2] = hv[n - 1][4] - hv[n - 2][1] * x[n - 1];
  505.     for (i = n - 3; i >= 0; i--)
  506.     x[i] = hv[i + 1][4] - hv[i][1] * x[i + 1] - hv[i][2] * x[i + 2];
  507.  
  508.     free(hv);
  509.     return TRUE;
  510. }
  511.  
  512.  
  513. /*
  514.  * Calculation of approximation cubic splines
  515.  * Input:  x[i], y[i], weights z[i]
  516.  *         
  517.  * Returns matrix of spline coefficients
  518.  */
  519. static spline_coeff *cp_approx_spline(plot, first_point, num_points)
  520. struct curve_points *plot;
  521. int first_point;        /* where to start in plot->points */
  522. int num_points;            /* to determine end in plot->points */
  523. {
  524.     spline_coeff *sc;
  525.     five_diag *m;
  526.     int xaxis = plot->x_axis;
  527.     int yaxis = plot->y_axis;
  528.     double *r, *x, *h, *xp, *yp;
  529.     /* HBB 980308: added 'GPHUGE' tag */
  530.     struct coordinate GPHUGE *this_points;
  531.     int i;
  532.  
  533.     sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff),
  534.                    "spline matrix");
  535.  
  536.     if (num_points < 4)
  537.     int_error("Can't calculate approximation splines, need at least 4 points", NO_CARET);
  538.  
  539.     this_points = (plot->points) + first_point;
  540.  
  541.     for (i = 0; i <= num_points - 1; i++)
  542.     if (this_points[i].z <= 0)
  543.         int_error("Can't calculate approximation splines, all weights have to be > 0", NO_CARET);
  544.  
  545.     m = (five_diag *) gp_alloc((num_points - 2) * sizeof(five_diag), "spline help matrix");
  546.  
  547.     r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
  548.     x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
  549.     h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
  550.  
  551.     xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
  552.     yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
  553.  
  554.     /* KB 981107: With logarithmic axis first convert back to linear scale */
  555.  
  556.     if (log_array[xaxis]) {
  557.     for (i = 0; i <= num_points - 1; i++)
  558.         xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
  559.     } else {
  560.     for (i = 0; i <= num_points - 1; i++)
  561.         xp[i] = this_points[i].x;
  562.     }
  563.     if (log_array[yaxis]) {
  564.     for (i = 0; i <= num_points - 1; i++)
  565.         yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
  566.     } else {
  567.     for (i = 0; i <= num_points - 1; i++)
  568.         yp[i] = this_points[i].y;
  569.     }
  570.  
  571.     for (i = 0; i <= num_points - 2; i++)
  572.     h[i] = xp[i + 1] - xp[i];
  573.  
  574.     /* set up the matrix and the vector */
  575.  
  576.     for (i = 0; i <= num_points - 3; i++) {
  577.     r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
  578.             - (yp[i + 1] - yp[i]) / h[i]);
  579.  
  580.     if (i < 2)
  581.         m[i][0] = 0;
  582.     else
  583.         m[i][0] = 6 / this_points[i].z / h[i - 1] / h[i];
  584.  
  585.     if (i < 1)
  586.         m[i][1] = 0;
  587.     else
  588.         m[i][1] = h[i] - 6 / this_points[i].z / h[i] * (1 / h[i - 1] + 1 / h[i])
  589.         - 6 / this_points[i + 1].z / h[i] * (1 / h[i] + 1 / h[i + 1]);
  590.  
  591.     m[i][2] = 2 * (h[i] + h[i + 1])
  592.         + 6 / this_points[i].z / h[i] / h[i]
  593.         + 6 / this_points[i + 1].z * (1 / h[i] + 1 / h[i + 1]) * (1 / h[i] + 1 / h[i + 1])
  594.         + 6 / this_points[i + 2].z / h[i + 1] / h[i + 1];
  595.  
  596.     if (i > num_points - 4)
  597.         m[i][3] = 0;
  598.     else
  599.         m[i][3] = h[i + 1] - 6 / this_points[i + 1].z / h[i + 1] * (1 / h[i] + 1 / h[i + 1])
  600.         - 6 / this_points[i + 2].z / h[i + 1] * (1 / h[i + 1] + 1 / h[i + 2]);
  601.  
  602.     if (i > num_points - 5)
  603.         m[i][4] = 0;
  604.     else
  605.         m[i][4] = 6 / this_points[i + 2].z / h[i + 1] / h[i + 2];
  606.     }
  607.  
  608.     /* solve the matrix */
  609.     if (!solve_five_diag(m, r, x, num_points - 2)) {
  610.     free(h);
  611.     free(x);
  612.     free(r);
  613.     free(m);
  614.     free(xp);
  615.     free(yp);
  616.     int_error("Can't calculate approximation splines", NO_CARET);
  617.     }
  618.     sc[0][2] = 0;
  619.     for (i = 1; i <= num_points - 2; i++)
  620.     sc[i][2] = x[i - 1];
  621.     sc[num_points - 1][2] = 0;
  622.  
  623.     sc[0][0] = yp[0] + 2 / this_points[0].z / h[0] * (sc[0][2] - sc[1][2]);
  624.     for (i = 1; i <= num_points - 2; i++)
  625.     sc[i][0] = yp[i] - 2 / this_points[i].z *
  626.         (sc[i - 1][2] / h[i - 1]
  627.          - sc[i][2] * (1 / h[i - 1] + 1 / h[i])
  628.          + sc[i + 1][2] / h[i]);
  629.     sc[num_points - 1][0] = yp[num_points - 1]
  630.     - 2 / this_points[num_points - 1].z / h[num_points - 2]
  631.     * (sc[num_points - 2][2] - sc[num_points - 1][2]);
  632.  
  633.     for (i = 0; i <= num_points - 2; i++) {
  634.     sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
  635.         - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
  636.     sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
  637.     }
  638.  
  639.     free(h);
  640.     free(x);
  641.     free(r);
  642.     free(m);
  643.     free(xp);
  644.     free(yp);
  645.  
  646.     return (sc);
  647. }
  648.  
  649.  
  650. /*
  651.  * Calculation of cubic splines
  652.  *
  653.  * This can be treated as a special case of approximation cubic splines, with
  654.  * all weights -> infinity.
  655.  *         
  656.  * Returns matrix of spline coefficients
  657.  */
  658. static spline_coeff *cp_tridiag(plot, first_point, num_points)
  659. struct curve_points *plot;
  660. int first_point, num_points;
  661. {
  662.     spline_coeff *sc;
  663.     tri_diag *m;
  664.     int xaxis = plot->x_axis;
  665.     int yaxis = plot->y_axis;
  666.     double *r, *x, *h, *xp, *yp;
  667.     /* HBB 980308: added 'GPHUGE' tag */
  668.     struct coordinate GPHUGE *this_points;
  669.     int i;
  670.  
  671.     if (num_points < 3)
  672.     int_error("Can't calculate splines, need at least 3 points", NO_CARET);
  673.  
  674.     this_points = (plot->points) + first_point;
  675.  
  676.     sc = (spline_coeff *) gp_alloc((num_points) * sizeof(spline_coeff), "spline matrix");
  677.     m = (tri_diag *) gp_alloc((num_points - 2) * sizeof(tri_diag), "spline help matrix");
  678.  
  679.     r = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline right side");
  680.     x = (double *) gp_alloc((num_points - 2) * sizeof(double), "spline solution vector");
  681.     h = (double *) gp_alloc((num_points - 1) * sizeof(double), "spline help vector");
  682.  
  683.     xp = (double *) gp_alloc((num_points) * sizeof(double), "x pos");
  684.     yp = (double *) gp_alloc((num_points) * sizeof(double), "y pos");
  685.  
  686.     /* KB 981107: With logarithmic axis first convert back to linear scale */
  687.  
  688.     if (log_array[xaxis]) {
  689.     for (i = 0; i <= num_points - 1; i++)
  690.         xp[i] = exp(this_points[i].x * log_base_array[xaxis]);
  691.     } else {
  692.     for (i = 0; i <= num_points - 1; i++)
  693.         xp[i] = this_points[i].x;
  694.     }
  695.     if (log_array[yaxis]) {
  696.     for (i = 0; i <= num_points - 1; i++)
  697.         yp[i] = exp(this_points[i].y * log_base_array[yaxis]);
  698.     } else {
  699.     for (i = 0; i <= num_points - 1; i++)
  700.         yp[i] = this_points[i].y;
  701.     }
  702.  
  703.     for (i = 0; i <= num_points - 2; i++)
  704.     h[i] = xp[i + 1] - xp[i];
  705.  
  706.     /* set up the matrix and the vector */
  707.  
  708.     for (i = 0; i <= num_points - 3; i++) {
  709.     r[i] = 3 * ((yp[i + 2] - yp[i + 1]) / h[i + 1]
  710.             - (yp[i + 1] - yp[i]) / h[i]);
  711.  
  712.     if (i < 1)
  713.         m[i][0] = 0;
  714.     else
  715.         m[i][0] = h[i];
  716.  
  717.     m[i][1] = 2 * (h[i] + h[i + 1]);
  718.  
  719.     if (i > num_points - 4)
  720.         m[i][2] = 0;
  721.     else
  722.         m[i][2] = h[i + 1];
  723.     }
  724.  
  725.     /* solve the matrix */
  726.     if (!solve_tri_diag(m, r, x, num_points - 2)) {
  727.     free(h);
  728.     free(x);
  729.     free(r);
  730.     free(m);
  731.     free(xp);
  732.     free(yp);
  733.     int_error("Can't calculate cubic splines", NO_CARET);
  734.     }
  735.     sc[0][2] = 0;
  736.     for (i = 1; i <= num_points - 2; i++)
  737.     sc[i][2] = x[i - 1];
  738.     sc[num_points - 1][2] = 0;
  739.  
  740.     for (i = 0; i <= num_points - 1; i++)
  741.     sc[i][0] = yp[i];
  742.  
  743.     for (i = 0; i <= num_points - 2; i++) {
  744.     sc[i][1] = (sc[i + 1][0] - sc[i][0]) / h[i]
  745.         - h[i] / 3 * (sc[i + 1][2] + 2 * sc[i][2]);
  746.     sc[i][3] = (sc[i + 1][2] - sc[i][2]) / 3 / h[i];
  747.     }
  748.  
  749.     free(h);
  750.     free(x);
  751.     free(r);
  752.     free(m);
  753.     free(xp);
  754.     free(yp);
  755.  
  756.     return (sc);
  757. }
  758.  
  759. static void do_cubic(plot, sc, first_point, num_points, dest)
  760. struct curve_points *plot;    /* still containes old plot->points */
  761. spline_coeff *sc;        /* generated by cp_tridiag */
  762. int first_point;        /* where to start in plot->points */
  763. int num_points;            /* to determine end in plot->points */
  764. struct coordinate *dest;    /* where to put the interpolated data */
  765. {
  766.     double xdiff, temp, x, y;
  767.     int i, l;
  768.     int xaxis = plot->x_axis;
  769.     int yaxis = plot->y_axis;
  770.     /* HBB 980308: added 'GPHUGE' tag */
  771.     struct coordinate GPHUGE *this_points;
  772.  
  773.     /* min and max in internal (eg logged) co-ordinates. We update
  774.      * these, then update the external extrema in user co-ordinates
  775.      * at the end.
  776.      */
  777.  
  778.     double ixmin, ixmax, iymin, iymax;
  779.     double sxmin, sxmax, symin, symax;    /* starting values of above */
  780.  
  781.     if (log_array[xaxis]) {
  782.     ixmin = sxmin = log(min_array[xaxis]) / log_base_array[xaxis];
  783.     ixmax = sxmax = log(max_array[xaxis]) / log_base_array[xaxis];
  784.     } else {
  785.     ixmin = sxmin = min_array[xaxis];
  786.     ixmax = sxmax = max_array[xaxis];
  787.     }
  788.  
  789.     if (log_array[yaxis]) {
  790.     iymin = symin = log(min_array[yaxis]) / log_base_array[yaxis];
  791.     iymax = symax = log(max_array[yaxis]) / log_base_array[yaxis];
  792.     } else {
  793.     iymin = symin = min_array[yaxis];
  794.     iymax = symax = max_array[yaxis];
  795.     }
  796.  
  797.  
  798.     this_points = (plot->points) + first_point;
  799.  
  800.     l = 0;
  801.  
  802.     xdiff = (this_points[num_points - 1].x - this_points[0].x) / (samples - 1);
  803.  
  804.     for (i = 0; i < samples; i++) {
  805.     x = this_points[0].x + i * xdiff;
  806.  
  807.     while ((x >= this_points[l + 1].x) && (l < num_points - 2))
  808.         l++;
  809.  
  810.     /* KB 981107: With logarithmic x axis the values were converted back to linear  */
  811.     /* scale before calculating the coefficients. Use exponential for log x values. */
  812.  
  813.     if (log_array[xaxis]) {
  814.         temp = exp(x * log_base_array[xaxis]) - exp(this_points[l].x * log_base_array[xaxis]);
  815.         y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
  816.     } else {
  817.         temp = x - this_points[l].x;
  818.         y = ((sc[l][3] * temp + sc[l][2]) * temp + sc[l][1]) * temp + sc[l][0];
  819.     }
  820.     /* With logarithmic y axis, we need to convert from linear to log scale now. */
  821.     if (log_array[yaxis]) {
  822.         if (y > 0.)
  823.         y = log(y) / log_base_array[yaxis];
  824.         else
  825.         y = symin - (symax - symin);
  826.     }
  827.     dest[i].type = INRANGE;
  828.     STORE_AND_FIXUP_RANGE(dest[i].x, x, dest[i].type, ixmin, ixmax, auto_array[xaxis], NOOP, continue);
  829.     STORE_AND_FIXUP_RANGE(dest[i].y, y, dest[i].type, iymin, iymax, auto_array[yaxis], NOOP, NOOP);
  830.  
  831.     dest[i].xlow = dest[i].xhigh = dest[i].x;
  832.     dest[i].ylow = dest[i].yhigh = dest[i].y;
  833.  
  834.     dest[i].z = -1;
  835.  
  836.     }
  837.  
  838.     UPDATE_RANGE(ixmax > sxmax, max_array[xaxis], ixmax, xaxis);
  839.     UPDATE_RANGE(ixmin < sxmin, min_array[xaxis], ixmin, xaxis);
  840.     UPDATE_RANGE(iymax > symax, max_array[yaxis], iymax, yaxis);
  841.     UPDATE_RANGE(iymin < symin, min_array[yaxis], iymin, yaxis);
  842.  
  843. }
  844.  
  845. /*
  846.  * This is the main entry point used. As stated in the header, it is fine,
  847.  * but I'm not too happy with it.
  848.  */
  849.  
  850. void gen_interp(plot)
  851. struct curve_points *plot;
  852. {
  853.  
  854.     spline_coeff *sc;
  855.     double *bc;
  856.     struct coordinate *new_points;
  857.     int i, curves;
  858.     int first_point, num_points;
  859.  
  860.     curves = num_curves(plot);
  861.     new_points = (struct coordinate *) gp_alloc((samples + 1) * curves * sizeof(struct coordinate), "interpolation table");
  862.  
  863.     first_point = 0;
  864.     for (i = 0; i < curves; i++) {
  865.     num_points = next_curve(plot, &first_point);
  866.     switch (plot->plot_smooth) {
  867.     case CSPLINES:
  868.         sc = cp_tridiag(plot, first_point, num_points);
  869.         do_cubic(plot, sc, first_point, num_points,
  870.              new_points + i * (samples + 1));
  871.         free(sc);
  872.         break;
  873.     case ACSPLINES:
  874.         sc = cp_approx_spline(plot, first_point, num_points);
  875.         do_cubic(plot, sc, first_point, num_points,
  876.              new_points + i * (samples + 1));
  877.         free(sc);
  878.         break;
  879.  
  880.     case BEZIER:
  881.     case SBEZIER:
  882.         bc = cp_binomial(num_points);
  883.         do_bezier(plot, bc, first_point, num_points,
  884.               new_points + i * (samples + 1));
  885.         free((char *) bc);
  886.         break;
  887.     default:        /* keep gcc -Wall quiet */
  888.         ;
  889.     }
  890.     new_points[(i + 1) * (samples + 1) - 1].type = UNDEFINED;
  891.     first_point += num_points;
  892.     }
  893.  
  894.     free(plot->points);
  895.     plot->points = new_points;
  896.     plot->p_max = curves * (samples + 1);
  897.     plot->p_count = plot->p_max - 1;
  898.  
  899.     return;
  900. }
  901.  
  902. /* 
  903.  * sort_points
  904.  *
  905.  * sort data succession for further evaluation by plot_splines, etc.
  906.  * This routine is mainly introduced for compilers *NOT* supporting the
  907.  * UNIX qsort() routine. You can then easily replace it by more convenient
  908.  * stuff for your compiler.
  909.  * (MGR 1992)
  910.  */
  911.  
  912. static int compare_points(p1, p2)
  913. struct coordinate *p1;
  914. struct coordinate *p2;
  915. {
  916.     if (p1->x > p2->x)
  917.     return (1);
  918.     if (p1->x < p2->x)
  919.     return (-1);
  920.     return (0);
  921. }
  922.  
  923. void sort_points(plot)
  924. struct curve_points *plot;
  925. {
  926.     int first_point, num_points;
  927.  
  928.     first_point = 0;
  929.     while ((num_points = next_curve(plot, &first_point)) > 0) {
  930.     /* Sort this set of points, does qsort handle 1 point correctly? */
  931.     qsort((char *) (plot->points + first_point), num_points,
  932.           sizeof(struct coordinate), (sortfunc) compare_points);
  933.     first_point += num_points;
  934.     }
  935.     return;
  936. }
  937.  
  938.  
  939. /*
  940.  * cp_implode() if averaging is selected this function computes the new
  941.  *              entries and shortens the whole thing to the necessary
  942.  *              size
  943.  * MGR Addendum
  944.  */
  945.  
  946. void cp_implode(cp)
  947. struct curve_points *cp;
  948. {
  949.     int first_point, num_points;
  950.     int i, j, k;
  951.     double x, y, sux, slx, suy, sly;
  952.     enum coord_type dot;
  953.  
  954.  
  955.     j = 0;
  956.     first_point = 0;
  957.     while ((num_points = next_curve(cp, &first_point)) > 0) {
  958.     k = 0;
  959.     for (i = first_point; i < first_point + num_points; i++) {
  960.         if (!k) {
  961.         x = cp->points[i].x;
  962.         y = cp->points[i].y;
  963.         sux = cp->points[i].xhigh;
  964.         slx = cp->points[i].xlow;
  965.         suy = cp->points[i].yhigh;
  966.         sly = cp->points[i].ylow;
  967.         dot = INRANGE;
  968.         if (cp->points[i].type != INRANGE)
  969.             dot = UNDEFINED;    /* This means somthing other than usual *//* just signal to check if INRANGE */
  970.         k = 1;
  971.         } else if (cp->points[i].x == x) {
  972.         y += cp->points[i].y;
  973.         sux += cp->points[i].xhigh;
  974.         slx += cp->points[i].xlow;
  975.         suy += cp->points[i].yhigh;
  976.         sly += cp->points[i].ylow;
  977.         if (cp->points[i].type != INRANGE)
  978.             dot = UNDEFINED;
  979.         k++;
  980.         } else {
  981.         cp->points[j].x = x;
  982.         cp->points[j].y = y / (double) k;
  983.         cp->points[j].xhigh = sux / (double) k;
  984.         cp->points[j].xlow = slx / (double) k;
  985.         cp->points[j].yhigh = suy / (double) k;
  986.         cp->points[j].ylow = sly / (double) k;
  987.         cp->points[j].type = INRANGE;
  988.         if (dot != INRANGE) {
  989.             if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
  990.             cp->points[j].type = OUTRANGE;
  991.             else if (!autoscale_ly) {
  992.             if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
  993.                 cp->points[j].type = OUTRANGE;
  994.             } else
  995.             cp->points[j].type = INRANGE;
  996.         }
  997.         j++;        /* next valid entry */
  998.         k = 0;        /* to read */
  999.         i--;        /* from this (-> last after for(;;)) entry */
  1000.         }
  1001.     }
  1002.     if (k) {
  1003.         cp->points[j].x = x;
  1004.         cp->points[j].y = y / (double) k;
  1005.         cp->points[j].xhigh = sux / (double) k;
  1006.         cp->points[j].xlow = slx / (double) k;
  1007.         cp->points[j].yhigh = suy / (double) k;
  1008.         cp->points[j].ylow = sly / (double) k;
  1009.         cp->points[j].type = INRANGE;
  1010.         if (dot != INRANGE) {
  1011.         if ((cp->points[j].x > xmax) || (cp->points[j].x < xmin))
  1012.             cp->points[j].type = OUTRANGE;
  1013.         else if (!autoscale_ly) {
  1014.             if ((cp->points[j].y > ymax) || (cp->points[j].y < ymin))
  1015.             cp->points[j].type = OUTRANGE;
  1016.         } else
  1017.             cp->points[j].type = INRANGE;
  1018.         }
  1019.         j++;        /* next valid entry */
  1020.     }
  1021.     /* insert invalid point to separate curves */
  1022.     if (j < cp->p_count) {
  1023.         cp->points[j].type = UNDEFINED;
  1024.         j++;
  1025.     }
  1026.     first_point += num_points;
  1027.     }                /* end while */
  1028.     cp->p_count = j;
  1029.     cp_extend(cp, j);
  1030. }
  1031.